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<N ■ Abstract 

> , 

^ ■ A nonperturbative lattice study of QCD at finite chemical potential is 

(y-v \ complicated due to the complex fermion determinant and the sign problem. 

Here we apply the method of stochastic quantization and complex Langevin 

f~^ ■ dynamics to this problem. We present results for U(l) and SU(3) one link 

C^ , models and QCD at finite chemical potential using the hopping expansion. 

The phase of the determinant is studied in detail. Even in the region 

where the sign problem is severe, we find excellent agreement between the 

Langevin results and exact expressions, if available. We give a partial 

rS ' understanding of this in terms of classical flow diagrams and eigenvalues 

C^ ■ of the Fokker-Planck equation. 
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1 Introduction 

A nonperturbative understanding of QCD at nonzero baryon density remains 
one of the outstanding problems in the theory of strong interactions. Besides the 
theoretical challenge, there is a clear phenomenological interest in pursuing these 
studies, due to the ongoing developments in heavy ion collision experiments, at 
RHIC, LHC and the planned FAIR facility at GSI. 

The standard nonperturbative tool to study quarks and gluons, lattice QCD, 
cannot be applied in a straightforward manner, because the complexity of the 
fermion determinant prohibits the use of approaches based on importance sam- 
pling. This is commonly referred to as the sign problem. Since the start of the 
millenium a number of new methods has been devised. These include reweight- 
ing [B [21 [3l H] , Taylor series expansion in /x/T [5], [6], [71 [8] , imaginary chemical 
potential and analytical continuation [9l [TOl [Til [12] , and the use of the canon- 
ical ensemble [HI HI] and the density of states [15]. Except for the last two, 
these approaches are approximate by construction and aimed at describing the 
QCD phase diagram in the region of small chemical potential and temperatures 



around the crossover between the confined and the deconfined phase. In this pa- 
per we discuss an approach which is manifestly independent of those hsted above: 
stochastic quantization and complex Langevin dynamics. How well this method 
will work is not known a priori. However, one of the findings of our study is that 
excellent agreement is found in the case of simple models, where comparison with 
results obtained differently is available. In particular we find that the range of 
applicability is not restricted to small chemical potential and, importantly, does 
not depend on the severity of the sign problem. The first results we present for 
lattice QCD at nonzero density are encouraging, although a systematic analysis 
has not yet been performed. 

This paper is organized as follows. In the following section we briefiy describe 
the idea behind stochastic quantization and the necessity to use complex Langevin 
dynamics in the case of nonzero chemical potential. In Sees. [3] and |4] we apply 
this technique to U(l) and SU(3) one link models. In both cases a comparison 
with exact results can be made. We study the phase of the determinant in detail. 
In the case of the U(l) model, we employ the possibility to analyse classical fiow 
diagrams and the Fokker-Planck equation to gain further understanding of the 
results. In Sec. [5] we turn to QCD, using the full gauge dynamics but treating 
the fermion determinant in the hopping expansion. Our findings and outlook to 
the future are summarized in Sec. El The Appendix contains a brief discussion 
of the Fokker-Planck equation in Minkowski time for the one link U(l) model. 

2 Stochastic quantization and complex Langevin 
dynamics 

The main idea of stochastic quantization [T6l [T7] is that expectation values are 

obtained as equilibrium values of a stochastic process. To implement this, the 

system evolves in a fictitious time direction 6', subject to stochastic noise, i.e. it 

evolves according to Langevin dynamics. Consider for the moment a real scalar 

field 0(x) in d dimensions with a real euclidean action S. The Langevin equation 

reads 

d(t){x,d) (55 [0] 



de 5(p{x,e) 

where the noise satisfies 



+ r/(x,^), (2.1^ 



{7]{x,e)) = o, {r]{x,e)7]{x',e')) = 26{x-x')6{e-e'). (2.2) 

By equating expectation values, defined as 

{0[4>{x,9)]),= [d4>P[(P,9]0[^{x)], (2.3) 



where O is an arbitrary operator and the brackets on the left-hand side denote a 
noise average, it can be shown that the probability distribution P[0, 6] satisfies 



the Fokker-Planck equation 



3^")- /-A^^f^^ + JM |P|^.,], (2,4) 



do J 6(f){x,e) \6(j){x,e) 6(j){x,e) 

In the case of a real action S, the stationary solution of the Fokker-Planck equa- 
tion, P[0] ~ exp (— S'[0]), will be reached in the large time limit 9 -^ oo, ensuring 
convergence of the Langevin dynamics to the correct equilibrium distribution. 
When the action is complex, as is the case in QCD at finite chemical potential, 
the situation is not so easy. It is still possible to consider Langevin dynamics 
based on Eq. (12.1 p [TSl [T^ [201 EI]- However, due to the complex force on the 
right-hand side, fields will now be complex as well: <^ — > Re<^-|-ilm0. As a result, 
proofs of the convergence towards the (now complex) weight e~^ are problematic. 

In the past, complex Langevin dynamics has been applied to effective three- 
dimensional spin models with complex actions, related to lattice QCD at finite yU 
in the limit of strong coupling and large fermion mass [22l [23l [21] (for applications 
to other models, see e.g. Ref. [25]). Our work has also partly been inspired by 
the recent application of stochastic quantization to solve nonequilibrium quan- 
tum field dynamics [261 123 [21] • In that case the situation is even more severe. 
Nevertheless, numerical convergence towards a solution can be obtained under 
certain conditions, both for simple models and four-dimensional field theories. 
For an illustration we present some original results in the appendix. 

Here we consider models with a partition function whose form is motivated by 
or derived from QCD at finite chemical potential. The QCD partition function 
reads 

Z= [ DUe-^'^detM, (2.5) 



where Sb{U) is the bosonic action depending on the gauge links U and det M is 
the complex fermion determinant, satisfying 

det M(/i) = [det M(-^)]*. (2.6) 

Specifically, for Wilson fermions the fermion matrix has the schematic form 
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i=l 



(2.7) 
Here T are lattice translations, r±^ = 1 ± 7^, and k, is the hopping parameter. 
Chemical potential is introduced by multiplying the temporal links in the forward 
(backward) direction with e'^ (e~^) [29]. We use Eq. (12.71) as a guidance to 
construct the U(l) and SU(3) one link models considered next. 



3 One link U(l) model 

3.1 Complex Langevin dynamics 

We consider a one link model with one degree of freedom, written a.s U = e 
The partition function is written suggestively as 



/n dr 
dUe-^^detM= / — e"^^ det M, (3.1) 

where the "bosonic" part of the action reads 

S^(x) = --{U + U-^) =-pcosx, (3.2) 

while the "determinant" is constructed by multiplying the forward (backward) 
link with e^ {e~^), 

detM =l + -K [e^'U + e"''?/"^] = 1 + k cos(x - ifi). (3.3) 

Due to the chemical potential, the determinant is complex and has the same 
property as the fermion determinant in QCD, i.e. detM(/i) = [det M(— /i)]*. For 
an imaginary chemical potential /x = i/ij, the determinant is real, as is the case 
in QCD. 

Observables are defined as 

1 n dr 
{0{x)) = ^ —e~^^detMO{x). (3.4) 

^ J-7T 27r 

In this model most expectation values can be evaluated analytically. We consider 
here the following observables: 

• Polyakov loop: 

(^7) = (e'") = ^ [h{(3) + kI[{(3) cosh/i - «:Ji(/5)//5sinh/i] , (3.5) 

where the partition function equals 

Z = Io{p) + Kh{p) cosh/i, (3.6) 

and In{P) are the modified Bessel functions of the first kind. 

• Conjugate Polyakov loop: 

(f/-i) = (e-^^) = (e^^) . (3.7) 

At finite chemical potential, (U) and {U~^) are both real, but different. 



Plaquette: 



d 1 

(cosx) = — InZ = — [Ji(/5) + kI[{[3) cosh/i] . 



(3i 



Note that (cosx) = ^(e"' + e "'). 



Density: 



d / msm(x-i/i) \ 1 . /^ • i 

^^^ = 9;;^"^ = \ l + .cos(x -./.) / = z"'^(^)^^"^^- 



(3.9) 



At small chemical potential (n) increases linearly with /x, while at large 
chemical potential {n) — * 1 exponentially fast. 

We now aim to estimate these observables using numerical techniques. Due 
to the complexity of the determinant, they cannot be estimated using methods 
based on importance sampling. Instead, we attempt to obtain expectation values 
using stochastic quantization. 

At nonzero chemical potential, the action is complex and it becomes necessary 
to consider complex Langevin dynamics. We write therefore x ^ z = x + iy, and 
consider the following complex Langevin equations 



Xn+l — Xn + eKxyXn, lln) + V^Vn, 
yn+l =yn + eKy{Xn, Z/n)- 

Here we have discretized Langevin time as 6* = ne, and the noise satisfies 



iVn) = 0, 

The drift terms are given by 

dS,s 



{VnVn') = 2(5„ 



K, 



-Re- 



dx 



where the effective action reads 






x— >a:+«y 



S'eflf = S'b — In det M = —[3 cos x — In [1 + k cos(x — ijj)] . 
Explicitly, the drift terms are 

cosh(y — /i) + K cos X 



Kr = — sin X 



/3coshy + K- 



Ky = —(3 COS X sinh y — k smh{y — /i) 



D{x) 

cos X + K cosh(y — /i) 



where 



D{x) 

2 , r • • 1 / m2 



D{x) = [I + Kcosxcosh.{y — j^)] + [k sin x sinh (y — //)] . 



(3.10) 
(3.11) 



(3.12) 

(3.13) 

(3.14) 

(3.15) 
(3.16) 

(3.17) 




Figure 1: Real part of the Polyakov loop (e*^) (left) and the conjugate Polyakov 
loop (e""') (right) as a function of fi for three values of /3 at fixed n = 1/2. The 
lines are the analytical results, the symbols are obtained with complex Langevin 
dynamics. 



Occasionally we will also consider this model by expanding in small n, the hopping 
expansion, and take 



(3.18) 



Scs = —(3 cos a; — kcos{x — ifi) (hopping expansion). 



This limit is motivated by the model of Heavy Dense Matter used in Ref. [SB]. A 
direct application of our method to QCD in the hopping expansion is presented 
in Sec. |5l 

In order to compute expectation values, also the observables have to be com- 
plexified. For example, after complexification x ^ z = x + iy, the plaquette 
reads 

cosx — > cos(x + iy) = cos x cosh y — i sin x sinh ?/. (3.19) 

All operators we consider are now complex, with the real (imaginary) part being 
even (odd) under x — > —x. 

The Langevin dynamics can be solved numerically. In Fig. [1] the real parts 
of the Polyakov loop and the conjugate Polyakov loop are shown as a function 
of /i for three values of f5 at fixed k = 1/2. In Fig. [2] (left) the density is shown. 
The lines are the exact analytical results. The symbols are obtained with the 
stochastic quantization. We observe excellent agreement between the analytical 
and numerical results. For the results shown here and below, we have used 
Langevin stepsize e = 5 x 10~^ and 5 x 10^ time steps. Errors are estimated 
with the jackknife procedure. The imaginary part of all observables shown here 
is consistent with zero within the error in the Langevin dynamics^ This can be 
understood from the symmetries of the drift terms and the complexified operators, 



^Analytically they are identically zero. 




Figure 2: Left: Real part of the density (n). Right: Real part of the plaquette 
(cosx) versus iJ?. Results at positive (negative) /i^ have been obtained with 
complex (real) Langevin evolution. 



since the drift terms behave under x 



-X as 



K^{-x, y; /i) = -K^{x, y; fx), Ky{-x, y; /i) = Ky{x, y; //), 



(3.20) 



while the imaginary parts are odd. Therefore, after averaging over the Langevin 
trajectory the expectation value is expected to reach zero within the error, which 
is what we observe. As an aside, we note that the symmetry of the drift terms 
under y ^ —y, 



Kx{x,-y; 



-jj) = Kx{x,y;n), Ky{x,-y]-n) = -Ky{x,y; n), (3.21) 



relates positive and negative chemical potential. 

At imaginary chemical potential /i = ifij the determinant is real, so that the 
complexification of the Langevin dynamics is not necessary. We demonstrate 
the smooth connection for results obtained at imaginary /i using real Langevin 
dynamics and results obtained at real n using complex Langevin dynamics for the 
expectation value of the plaquette (cosx) in Fig. [2] (right). Since the plaquette 
is even under /i —^ — /i, we show the result as a function of /x^, so that the 
left side of the plot corresponds to imaginary chemical potential, while the right 
side corresponds to real chemical potential. On both sides excellent agreement 
with the analytical expression can be observed. We also note that the errors are 
comparable on both sides. 



3.2 Phase of the determinant 

At finite chemical potential the determinant is complex and can be written as 

detM= IdetMle*"^. (3.22) 
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-/i)). Right: Real part 



In order to assess the severity of the sign problem, we consider the phase of the 
determinant and study the behaviour of e**^. An observable often used for this 
purpose [381 EH] is 

det M{fi) \ I det M{ii) 



{e'"f) 



(3.23) 



, det M(/i)*/ \ det M(-/i) 

where we used Eq. f l2.6p . At zero chemical potential the ratio equals one, while 
at large /i one finds in this model that 

HP) 



lim (e 

/^— >oo 



{e^'^) 



im 



+ 0{e 



(3.24) 



for nonzero (3. In expressing Eq. (I3.23P as the expectation value obtained from 
the complex Langevin process, complex conjugation has to be performed after the 
complexification of the variables, as discussed above. In that case det M(/i)* as 
a complex number is not the complex conjugate of det M{fi). To avoid confusion 
we write det M(— /x) in all relevant expressions. Notice that this implies that (j) 
itself is also complex. 

In Fig. [3] (left) we show the real part of this observable as a function of /x. The 
imaginary part is again zero analytically and zero within the error in the Langevin 
process. The lines are obtained by numerical integration of the observable with 
the complex weight, while the symbols are again obtained from Langevin dy- 
namics. We note again excellent agreement between the semi-analytical and the 
stochastic results. In particular, there seems to be no problem in accessing the 
region with larger fi where the average phase factor becomes very small. The 
numerical error is under control in the entire range. We find therefore that the 
sign problem does not appear to be a problem for this method in this modelo 



^In QCD, the average phase factor is expected to go to zero exponentially fast in the ther- 
modynamic limit. 
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Figure 4: Scatter plot of e^*"^ = detM(/i)/ det M(—/i) during the Langevin evo- 
lution for various values of /x at /9 = 1, k = 1/2. Note the different scale in the 
middle box. 
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In contrast to what could be inferred from the result above, expectation values 
of e^'^ are not phase factors, due to the complexity of the action. This can be 
seen by considering 



_2,;^ /detM(-ri\ _Z(-fi) 



where the second and third equality follow from the cancelation of det M{fi) 
in the definition of the expectation value and from Z being even in /i. We 
have also computed this observable using Langevin dynamics and the result is 
shown in Fig. [3] (right). For the Langevin parameters used here, we observe that 
the numerical estimate is consistent with 1, but with quite large errors when 
fi increases at small values of (3. We have found that increasing the Langevin 
time reduces the uncertainty. We conclude that at large chemical potential this 
ratio of determinants is the most sensitive and slowest converging observable we 
considered. 

We observed above that the average phase factor becomes very small at large fi 
but that this does not manifest itself in all but one observable we consider. Insight 
into this feature can be gained by studying scatter plots of the phase factor during 
the Langevin process. In Fig.lHwe show the behaviour of e^**^ during the Langevin 
evolution in the two-dimensional plane spanned by Ree^*"^ and Ime^*'^. At zero 
chemical potential, Ree^**^ = 1 and Ime^*'^ = during the entire evolution. For 
increasing /i one finds more and more deviations from this, with an interesting 
structure at intermediate values of /i. Note that the resulting distribution is 
approximately invariant under reflection in Ime^*'^ — > — Ime^*'^, ensuring that the 
imaginary part of the expectation value (e^**^) vanishes within the error. Due to 
the wide distribution, the horizontal and vertical scales in the middle section of 
Fig. m are much larger than in the top and bottom part. However, the average 
phase factor remains well defined for all values of /i, as can be seen in Fig. [31 At 
large /i, the average phase factor becomes very small. However, the distribution 
is very narrow, see Fig. H] (bottom). Therefore, although the average is close to 
zero, the error in the Langevin dynamics is very well under control. 

3.3 Fixed points and classical flow 

The excellent results obtained above can partly be motivated by the structure of 
the dynamics in the classical limit, i.e. in absence of the noise. As we demonstrate 
below, the classical fiow and fixed point structure is easy to understand when 
fi = and, most importantly, does not change qualitatively in the presence of 
nonzero chemical potential. 

Classical fixed points are determined by the extrema of the classical (effective) 
action, i.e. by putting K^ = Ky = 0. We first consider the "bosonic" model and 
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take K = 0. The drift terms arqj 

Kx{x,y) = —P sin X cosh y, Ky{x, y) = —P cosxsmhy. (3.26) 

We see that there is one stable fixed point at (x, y) = (0, 0) and one unstable 
fixed point at (vr, 0). Moreover, the classical fiow equation, dz/dO = —(3smz, can 
be solved analytically, with the solution 

tan^ = e-^(''"^«)tan^, (3.27) 

2 2 ' ^ ' 

where ^(0) is the initial value at ^ = ^o- We find therefore that the stable fixed 
point is reached for all z{0), except when x{0) = tt. On this line the solution 
reads 

tanh ^ = e'^^^-^") tanh ^, (3.28) 

2 2 ' ^ ' 

and the fiow diverges to y — > ±oo, except when starting precisely on the unstable 
fixed point (vr, 0). Note, however, that the noise in the x direction will kick the 
dynamics of the unstable trajectories. 

We now include the determinant, starting with the hopping expansion fl3.18p . 
Putting K^ = Ky = yields again one stable fixed point at {x,y) = {0,y^:) and 
one unstable fixed point at (tt,?/*), where 

Ksinhu 

tanhy, = — ^. 3.29 

p + K cosh /i 

Note that in the strong coupling limit y* = fi. We find therefore a simple modifi- 
cation of the bosonic model: in response to the chemical potential the two fixed 
points move in the vertical y direction, but not in the x direction. 

We continue with the full determinant included. Consider the case with /i = 
first, where real dynamics can be considered. Again we find the stable fixed point 
at X = and the unstable fixed point at x = tt. Provided that 

..(i + i)>l, (3.30) 

there are no additional fixed points. In order to satisfy this condition, we take 
K < 1 throughout. Using complex dynamics, while keeping /i = 0, we find that 
the stable fixed point at (x, y) = (0, 0) remains, but that there are now three 
unstable fixed points at x = vr, given by (vr, 0) and (vr, iy^,), with coshy* = 7. 
Interestingly, the fixed-point structure is therefore different for real and complex 
flow. 

Finally, we come to the full determinant at finite chemical potential. In this 
case the fixed points can only be determined numerically. We find a stable fixed 



•^For the bosonic model, there is of course no need to complexify the Langevin dynamics and 
one may take y = 0. This yields the same fixed points. 
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Figure 5: Classical flow diagram in the x — y plane for (3 = 1, /t= 1/2, /i = 0.1 
(top) and /i = 1 (bottom). The big dots indicate the fixed points at a; = and tt. 
The small circles indicate a trajectory during the Langevin evolution, each dot 
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Figure 6: As in the previous figure, with /i = 2 (top) and /i = 5 (bottom). 
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point at {x,y) = {0,ys) and unstable fixed points at {x,y) = (vr, ?/„). The y 
coordinates of tliese fixed points are determined by 



Ky 



Ksinhfy — u) , ^ 

= tP sinh y T t- TT^ = 0' 3.31 

x=0,TT 1 ± K COSIl(y — j2) 



where the upper (lower) sign applies to ys (yu)- At x = there is only one 
solution, while at a; = vr we find numerically that there are three (unstable) 
solutions for small chemical potential, two for intermediate /i and only one for 
large fi. Although the number of fixed points at x = vr depends on fi, we find 
that they are always unstable such that this has no effect on the dynamics, which 
is attracted to the stable fixed point at x = 0. 

In Figs. [5] and [6] we show the classical flow diagrams in the x — y plane. The 
direction of the arrows indicates {K^, Ky), evaluated at (x, y). The lengths of the 
arrows are normalized for clarity. The fixed points are indicated with the larger 
black dots. In the bosonic model (k = 0), the analytical solution showed that 
the fixed point at x = is globally attractive, except when x = vr. At nonzero k 
and fi, the fixed point at x = appears to be globally attractive as well, except 
again for x = vr. The small (blue) dots are part of a Langevin trajectory, each 
dot separated from the previous one by 500 steps. For vanishing fi, the dynamics 
takes place in the x direction only; for increasing /i it spreads more and more in 
the y direction. An interesting asymmetry around the classical fixed point in the 
y direction can be observed. However, the dynamics remains well contained in 
the X — y plane. 

We conclude therefore that the complex Langevin dynamics does not change 
qualitatively in the presence of a chemical potential, small or large. We take this 
as a strong indication that the method is insensitive to the sign problem. 

3.4 Fokker- Planck equation 

The microscopic dynamics of the Langevin equation. 



dx dS 

89 dx 



V, (3.32) 



where 6 is the (continuous) Langevin time, can be translated into the dynamics 
of a distribution P(x, 6), via the relation 

{0{x,9)), = J^P{x,9)0{x). (3.33) 

From the Langevin equation, it follows that P(x, 9) satisfies a Fokker-Planck 
equation, 

^P{x,9) = L'^pP{x,9), (3.34) 
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where Lpp is the complex Fokker-Planck operator 

L^ =^f^ + ^V (3.35) 

^^ dx \dx dx J 

The stationary solution of the Fokker-Planck equation is easily found by putting 
LppP(x) = and is given by Pst{x) ~ exp[— S'(x)]. In order to cast Eq. (13.341) 
as an eigenvalue problem, we write P{x,9) = e~^^P'^^\x). The solution of the 
Fokker-Planck equation can then be written as 

P(x,e) = ^— - + Ve-^V^)(a;)- (3.36) 

X 

It is therefore interesting to study the properties of the Fokker-Planck equation 
and the nonzero eigenvalues A. 

In order to do so, we consider the model in the hopping expansion (I3.18p . 
with the action 

S = —P cosx — Kcos{x — ifi). (3.37) 

Explicitly, the Fokker-Planck equation then reads 

P(x, 6) = P"{x, 9) + [/? sin X + K sin(x - i/i)] P'(x, 9) 

+ [[3cosx + Kcos(x — iji)] P(x, 9), (3.38) 

where primes/dots indicate x/9 derivatives. Using periodicity, P(x + 27r, 9) = 
P{x,9), we decompose 

— e'--p{x,9), (3.39) 

and we find 

P„(e) = -n^Pn{9) - nc+Pn+i{9) + nc_P„_i(e), (3.40) 

with 

c± = l{P + «:e±^) . (3.41) 

We note that this equation is completely real, such that all P„(^)'s are real. This 
is expected for the stationary solution, since from S*{x) = S{—x) it follows that 
Pg*t(x) = Pst{—x) and therefore P*^^ = P„,st- The numerical solution of Eq. (13.401) 
is shown in Fig. [7] for a number of modes Pn{9) for /i = 1 (left) and 3 (right). 
The initial values Pn(0) = 1 for all n. The number of modes is truncated, with 
—50 < n < 50. For large ±n, Pn{9) —>■ exponentially fast. The zero mode Pq 
is 9 independent and equal to 1 by construction. We have verified that the other 
modes converge to the values determined by the stationary solution ~ e""^. 
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Figure 7: Solution of the complex Fokker-Planck equation: Langevin time de- 
pendence of the modes P„(6') for various values of n, with /3 = 1, k = 1/2, and 
/i = 1 (left) and /i = 3 (right). 

The convergence properties can be understood from the eigenvalues of the 
Fokker-Planck operator. Writing Pn{0) = exp(— A^)P„ gives the eigenvalue equa- 
tion 

n'Pn + nc+Pn+i - nc_Pr,.i = XPn. (3.42) 

Since all Pn are real, also all eigenvalues A are real. Although at first sight 
this may seem surprising, it is a consequence of the symmetry of the action and 
therefore it also holds in e.g. the full model. 

The case A = corresponds to the stationary solution. Here we consider 
A 7^ 0. First take n = 0. It follows from Eq. (13^421) that Pq = 0. As a result, the 
sequences for n > and n < split in two. Written in matrix form, they read 



(3.43) 



V : : : : 

Approximating this matrix by a large but finite matrix, one can easily compute 
the eigenvalues numerically. We find that they are all positive and that the n > 0, 
n < sequences yield identical eigenvalues. In Fig. [H] (left) the four smallest 
nonzero eigenvalues are shown as a function of chemical potential. All eigenvalues 
are strictly positive and increase with /i. In Fig. [8] (right) the dependence on 
/5 is indicated. At vanishing (3, the /i dependence cancels, since in that case 
c+c_ = fi;^/4. Also as a function of /5 we observe that the eigenvalues are strictly 
positive. 

If the action and therefore the Langevin dynamics would be real, these results 
would be sufficient to explain the convergence of the observables towards the 
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Figure 8: Left: Four smallest nonzero eigenvalues of the complex Fokker-Planck 
equation as a function of n with (3 = 1, k = 1/2. Right: Smallest nonzero 
eigenvalue as a function of j3 for various values of /x at k = 1/2. 



correct values as observed above, by employing Eq. (13.361) in the large 6 limit. In 
the complex case we consider here, this is not immediately clear. Given the real 
Langevin equations. 



dx 

m 



Kx + r/. 



de "' 



(3.44) 



one can also consider the real distribution p{x, y, 6), satisfying the Fokker-Planck 

equation 

d 

-p{x,y,9) = LFPp{x,y,9), 



09' 
with the real Fokker-Planck operator 

d f d 



Lfp = t^ ( 7^ Kx] - T-^y 

ox \ox J oy 

After complexification, expectation values should then satisfy 

dxdy 



(3.45) 



(3.46) 



{Oix + ty,9)y^ 



2tt 



p{x,y,9)0{x + iy). 



(3.47) 



In contrast to the complex distribution P{x,9), the distribution p{x,y,9) is real 
and has the interpretation as a probability distribution in the x — y plane. The 
real and complex Fokker-Planck operators can be related, using 



— (0(x + ly, 9)) = f ^ 0(x + iy)LFPp{x, y, 9) 
o9 J 27r 

dxdy 
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■0{x^%y)L%^p{x,y,9). 



(3.48) 



Here partial integration for finite 9 is used as well as the analytic dependence 
oi O on z = X + iy. Eq. (13.481) suggests a relation between the spectrum of the 
complex and the real Fokker-Planck operator. However, we have not yet been 
able to show that a stationary solution of the real Fokker-Planck equation exists. 
We hope to come back to this issue in the future. 



4 One link SU(3) model 



4.1 Model 

In this section we consider a one link model where the link U is an element of 
SU(3). The partition function reads 

t/t/e-^^detM, (4.1) 

with the bosonic part of the actioiu 

^5 = -^(Trf/ + Trt/-^). (4.2) 

6 

For the fermion matrix we take 

M= l + fi:(e^a+[/ + e^'V_f/"^), (4.3) 

with (j± = |(1 ±(73). We use the Pauli matrix a^ rather than 7 matrices to avoid 
factors of 2. The determinant has the product form 

det M = det (1 + ne^'U) det (l + Ke-^f/"^) , (4.4) 

where the remaining determinants on the right-hand side are in colour space. In 
order to exponentiate the determinant, we use the identity, valid for U G SL(3, C), 

det {l + cU) = l + cTr U + c^Tr U~^ + cl (4.5) 

We find therefore that 

detM = e-^^, SF = -\nM^'''> -\nM^^\ (4.6) 

with the quark and anti-quark contributions 

yW(9) = 1 + Si^e^'P + ?,K^e^^p-^ + K^e^", (4.7) 

yW® = 1 + 3Ke-^p-i + 3K2e-2'^P + K^g-s^. (4.8) 



*Note that for an SU(3) matrix, U "^ ~ U^ . Nevertheless, we write U ^ to allow for a 
straightforward complexification of the Langevin dynamics. 
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Here we introduced the Polyakov loop P and its "conjugate" P ^, 



1. 



1. 



P = -Tr U, P~^ = -Tr U~\ (4.9) 

o o 

Note that PP^^ 7^ 1. At large /i, the anti-quark contribution A^^^^ -^ 1 and no 
longer contributes. However, the term is crucial to preserve the symmetry fl2.6p 
and is in particular relevant at imaginary and small real /i. 
Observables are defined as 

(0(f/)) = y f dU e"^«(^) det M{U) 0{U). (4.10) 

The observables we consider are the Polyakov loop P, the conjugate Polyakov 
loop P~^ and the density n. The latter is determined by 

d\nZ 



in) 



dn 



(4.11) 



and reads 



n 



ainX^") 91nA1® 



9/i 



dfi 



Kei'P + 2K^e^^'P-^ + K^e^^" Ke~^'P~^ + 2K'^e-^f'P + n^e-^'' , 
= 3 Ai(S 3 ^;jg^ . (4.12) 

At zero chemical potential, the density vanishes while at large [x the density 
77, — > 3, the maximal numbers of (spinless) quarks on the site. 

In this model, expectation values can be obtained directly by numerical inte- 
gration, allowing for a comparison with the results from stochastic quantization 
presented below. Since we only consider observables that depend on the conju- 
gacy class of f/, we only have to integrate over the conjugacy classes. These are 
parametrized by two angles — vr < 0i, 02 < tt. The reduced Haar measure on the 
conjugacy classes [f/] reads 



d{\J\ = ^ sin^ 



2^ 



sm 



i(^, + 



sm 



>= 



where A/" is a normalization constant. The matrix is parametrized as 

[/ = diag (e^-^S e^'^^ e-*('^i+'^^)) , 



such that 



S 



B 



- [cos(0i) + COS(02) + COS(01 + ^2)] 



and 
P 



[e' 



+ e 



-i{4>i+4'2) 



]. ^-' 



[^- 



J(<j}l+<j}2) 



(4.13) 

(4.14) 
(4.15) 

]. (4.16) 



It is now straightforward to compute expectation values by numerical integration 
over 01 and 02. 
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4.2 Complex Langevin dynamics 

In contrast to in the U(l) model, the Langevin dynamics is now defined in terms 
of matrix multiphcation. We denote U{6 + e) = U' and U{6) = U, where 6 is 
again the Langevin time and consider the Langevin process, 

U' = RU, R = exp [i\a [eKa + V^Va)] ■ (4.17) 

Here A^ (a = 1, . . . , 8) are the traceless, hermitian Gell-Mann matrices, normal- 
ized as Tr \a\b = 25ab. The noise satisfies 

{Tla) = 0, iVaVb) = 2Sab, (4.18) 

and the drift term reads 

Ka = -DaS,s, S,s = Sb + Sf. (4.19) 

Differentiation is defined as 

DafiU) = ^f (e^^^'-f/) . (4.20) 

In particular, DaU = iXaU and DaU~^ = ~iU~^Xa- 
The explicit expressions for the drift terms are 

K, = K^ + iff, (4.21) 

with 

Kf = -DaSB{U) = ^ [DaP + DaP-') , (4.22) 

K, = -D,Sf(U) = 3 j^^ + 3 j^^ , 

(4.23) 

written in terms of 

DaP = '-Tt XM, DaP-' = -^Tr fJ-^A,. (4.24) 

Let us first consider the case without chemical potential and take U G SU(3). 
Then it is easy to see that K^ = Ka and therefore K^R = 1. Moreover, since 
the Gell-man matrices are traceless, deti? = 1. Therefore, if U is an element of 
SU(3), it will remain in SU(3) during the Langevin process. The same results are 
found at finite imaginary chemical potential /i = ifij. At nonzero real n on the 
other hand, we find that R^R ^ 1, although deti? = 1 still holds. Therefore U 
will be an element of SL(3, C) during the Langevin evolution. If f/ is parametrized 
as ?7 = exp (iAaAa/2), this implies that the gauge potentials Aa are complex. 
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Figure 9: Real part of the Polyakov loop (P) (left) and the conjugate Polyakov 
loop {P~^) (right) as a function of fi for three values of (3 at fixed k, = 1/2. 
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Figure 10: Scatter plot of the Polyakov loop for j3 = 1, n = 1/2 and /i = 0, 1, 2, 3. 
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Figure 11: Left: Real part of the density (n). Right: Deviation from SU(3) 
during complex Langevin evolution: Tr U'^U /?> as a function of Langevin step, for 
/i= l,2,3,4at /3= l,ft= 1/2. 



We solved the Langevin process fl4.17p numerically. The matrix R is com- 
puted by exponentiating the complex traceless matrix i\a {eKa + y/eria), employ- 
ing Cardano's method [37] for finding the eigenvalues. In Fig. [9] the real part 
of the Polyakov loop (P) and the conjugate Polyakov loop (-P"^) are shown as 
a function of /i for three values of (3 at fixed k = 1/2. The lines are the 'ex- 
act' results obtained by numerically integrating over the angles (pi and 02, as 
discussed above. The symbols are obtained with complex Langevin dynamics, 
using the same Langevin stepsize and number of time steps as in the U(l) model 
(e = 5 X 10~^ and 5 x 10^ time steps). Errors are estimated with the jackknife 
procedure. Again, the imaginary part is zero analytically and consistent with 
zero within the error in the Langevin dynamics. Excellent agreement between 
the exact and the stochastic results can be seen. 

Scatter plots of the Polyakov loop during the Langevin evolution are shown in 
Fig. [To] for four values of // at /3 = 1 and k = 1/2. Every point is separated from 
the previous one by 500 time steps. Note that the distribution is approximately 
symmetric under reflection ImP — > — ImP, ensuring that Im (P) = within the 
error. We observe that the characteristic shape visible at // = becomes more 
and more fuzzy at larger fi, but the average remains well defined for all values 
of /i. The density is shown in Fig. [11] (left), with again good agreement between 
the exact and the stochastic results. We observe that saturation effects set in for 
smaller values of /i compared to the U(l) model, e.g. n = nmax/2 = 3/2 already 
at /i ?^ 1. 

During the complex Langevin evolution the matrix U G SL(3, C). In order to 
quantify how much it deviates from SU(3), we may follow the evolution of 



f{U) = -Ti U^U. 



(4.25) 
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Figure 12: Left: Real part of {e^'''') = (det M(/i)/det Af(-/i)). Right: Real part 
of (detM(-/i)/detM(/i)). 

It is easy to show that f{U) > 1, with the equality in the case that U E SU(3)|j 
It provides therefore a good measure to quantify the deviation from SU(3). In 
Fig. [11] (right), we show this quantity during the Langevin evolution. We observe 
that the deviations from 1 are present but not too large. If U is parametrized as 
U = exp {i\aAa/2), this implies the imaginary parts of the gauge potentials Aa 
do not become unbounded. 



4.3 Phase of the determinant 

As in the U(l) model, we study the phase of the determinant in the form 

det M{ij) 



{e'"^) 



det M{-fi) 



(4.26) 



At zero chemical potential, the ratio is 1. Due to the SU(3) structure, however, 
the behaviour at large /i is qualitatively different. We find 



lim det M{fi) = K^e^^ [l + Se"^ {k + k~^) P~^ + 0{t 



-2/^ 



/i-^oo 



lim det M{-fi) = n^e^'' [l + 3e~^ {k + k~^) P + 0{e 



2/^ 



/i— >oo 



such that 



det Mill) 
lim - — -— 

M^oo det iw (— /i) 



1 + Se^'^ [k + K^^) (P^i -P)+ 0{e 



2^i^ 



(4.27) 
(4.28) 

(4.29) 



^Consider U e SL(A^,C) and /([/) = Tr[/tt//7V. Using a polar decomposition, U = VP, 
with V S SU(A^) and P a positive semidcfinite hemiitian matrix with detP = 1, it is easy to 
show that /([/) > 1, with the equal sign holding when U G SU(A^). 
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Figure 13: Scatter plot of e^^'^ = det M(yu)/det M(— /x) during the Langevin 
evolution, for (3 = 1, k, = 1/2 and 0.25 < /i < 4. 
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As a result the average phase goes to 1 at large /i and not towards as in the U(l) 
modeljj Therefore we expect the sign problem to become exponentially small in 
the saturation regime at large /i. 

We have computed the average phase factor and the results are shown in 
Fig. [12] (left). The lines are again the 'exact' results. As is clear from this plot, 
the sign problem is quite mild for all values of /i, since the maximal deviation 
from 1 is less than 15%. In Fig. [H (right) the ratio (det M(-/i)/det M(/i)) = 
Z[—n)/Z{fi) is shown. Here we observe a small but systematic deviation from 
1, more pronounced at smaller f3 and intermediate /i. However, we found that 
the deviation from 1 is reduced when continuing the Langevin evolution to larger 
and larger times. As in the U(l) model, this observable is the most sensitive and 
slowest converging quantity. 

Scatter plots of the phase are presented in Fig. [131 At small chemical potential 
(top figure), there appears to be a similar structure as in the U(l) model, although 
not as pronounced. In the intermediate region (middle), the distribution is wider. 
At large fi (bottom), the distribution becomes narrow again, centered around 
(1,0). 

5 QCD at finite chemical potential 

5.1 Hopping expansion 

In this section we leave the one link models behind and consider QCD at chemical 
potential. The SU(3) gauge field contribution to the euclidean lattice action iqj 



Sb[U] = -/5 E E (^ tTr [/.,,. + Tr U'l^] - l] 
with (3 = G/g"^. The plaquettes are defined as 



(5.1) 



Ux,fiu — Ux,^Ux+fi^uU^j^f^^^U^j^, (5.2) 

and 

The fermion matrix M for Wilson fermions was already given in Eq. (12. 7p . The 
7 matrices satisfy 7^ = 7^ and 7^ = 1. We use periodic boundary conditions 
in space and antiperiodic boundary conditions in the euclidean time direction; 
the temperature and the number of time slices A^,- are related as T = l/N^ (the 
lattice spacing a = 1). The fermion matrix obeys 

Mt(/i) = 75M(-/x)75, (5.4) 



^This difference can be traced back to Eq 
^We write C/~^ rather than U'^] see footnote (4] 
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such that the determinant obeys Eq. (12. 6p . 

We consider the heavy quark hmit, where all spatial hopping terms are ignored 
and only the temporal links in the fermion determinant are preserved. We write 
therefore 

det M ^det[l- K {eT+^U^^^T^ + e-T.^U-lT.^)] 

= det (1 - 2Ke^U^,4T4)^ det (l - 2fi;e-''f/-]T„4)^ 

= JJ det (1 + /le^/^Px) ^ det (l + he-^^'^VZ^f , (5.5) 

X 

where we defined h = {2k)^^ and the (conjugate) Polyakov loops are 

A^r-l 

^x = n f^(^.-)-4' ^x ^ = n f^r.,x),4- (5-6) 

r = T=Nr~l 

In the final line of Eq. (15.51) the determinant refers to colour space only. The + 
sign appears because of the antiperiodic boundary conditions. 

This approximation is motivated by the Heavy Dense Model considered e.g. 
in Refs. [il[3ll[32l[33l[34l[35l[36], in which the limit 

K ^ 0, /U ^ oo, ne^ fixed, (5.7) 

was taken. However, here also the backward propagating quark, with the inverse 
Polyakov loop, is kept in order to preserve the relation (12. 6p . 
Using Eq. (14. 5p . the determinant can now be written as 

detM = e-^^, SF = -^{2\nM^^'^ + 2\nMf), (5.8) 

X 

with the quark and anti-quark contributions 

>l(^5) = 1 + 3hd'^^P^ + Sh^e^f'/^P^^ + h^e^f"/^, (5.9) 

7W(5) = 1 + 3he-^/^p-^ + S/i^e-^^/^Px + /i^e"^^/^, (5.10) 

where 

Px = ^TrPx, P-' = ^TtV:\ (5.11) 



The density is given by 



(iV) = E(-x)=T^, (5.12) 



27 



and we find 



n^ = 2T — ^^ + 2T 



6 






M 



(9) 



—6 7- . (5.13) 

At zero chemical potential, the density vanishes while at large /i the density 
rix -^ 6, the maximal numbers of quarks on a site. 

5.2 Complex Langevin dynamics 

The implementation of the Langevin dynamics follows closely the one discussed 
in the previous section on the SU(3) one link model. We denote Ux,^i{0 + e) = U'^^^ 
and Ux,^i{0) = Ux,^j., and consider the process 



^ x,^i J^x,^i 'J X,jJ.1 


^x,fi - 


= exp [iXa [eKx^^a + \^Vx^la)] , 


(5.14) 


with the noise satisfying 








(^x/.a) = 


= 0, 


iVxtiaVyub) = '^5^i,5ab5xy 


(5.15) 


The drift term is 










J^xfia ~ 


- -Dx,aS[U]. 


(5.16) 



Differentiation is defined as 

J-(TI\ = 

da 



Dx,af{U) = -^f (e^°'»f/,,^) . (5.17) 



a=0 



The drift term is written as 

J^x/ia = J^x/ia + ^xfiaJ (5.18) 

with the bosonic contribution 

Ka = -Dx,aSB[U] 

= ^^ 5^ Tr [XaUx,Cx,,u - Cx,^.MxiXa] , (5.19) 



where 



(^x,^J.u — Ux+fi,uU^^i^^^U^j^ + U,^_^_i^_~^^U^_p^^Ux-i>,u, (5.20) 

^x,fip = fJ x,u'J x+u,nf^ x+fi,i^ ' ^x-j>,j/^a;-!>,^tyx+/i-£',i" (5.21) 
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The ferinionic contribution is 



K'a = -D.,aSF[U] = 6,,Kl, (5.22) 



with 



K^„ = 6 r-, h 0- 



(5.23) 



The derivatives are 

T-l Nr-1 

D,aP. = D(r^)aP. = ^Tr J] ?7(,'x)4Aa n ^(-"-)4' (5-24) 

t'=0 t"=t 

T 

l^.aP-^ = /^(.x)aPx ' = -^Tr n ^(;'x)4Aa H V'x)4- (5-25) 

T<=N^-\ t"=t-1 

We have solved Eq. (I5.14p numerically. A detailed analysis is postponed to a 
future publication; here we present some results for illustration purposes. We 
have used the temporal gauge, where only the last link differs from the identity, 

f/(7V.-l,x)4 = Kc, f/(.,x)4 = l (Ty^N^-l). (5.26) 

To simplify the exponentiation we use the following updating factor in Eq. (I5.14p . 

Rx,^= n e*^''('^^'^"+^''^''''), (5.27) 

a=Pcrm(l,..,8) 

with random ordering from sweep to sweep and where K^^a is complex. R and 
R only differ by terms of order e^, which is the general systematic error of the 
Langevin algorithm. For the results shown here, we have employed Langevin 
stepsize e = 2 x 10~^ and 50000 iterations of 50 sweeps each, using ergodicity 
to calculate averages. Runaway trajectories have practically been eliminated 
by monitoring the drift and using adaptive step size. The lattice has size 4^, 
with /? = 5.6, K = 0.12. We have studied chemical potentials in the range 
0.5 < /i < 0.9, using A'^^ = 3 fermion flavours. 

In Fig. [13] we present the real part of the Polyakov loop and the conjugate 
Polyakov loop (left) and the density (right). These results appear consistent with 
those obtained in Ref. [36] using reweighting techniques, although at this level of 
the study both statistics and thermalization are not yet optimal. Nevertheless 
we clearly see that at /i = 0.5 the system is in the low-density "confining" phase 
whereas for larger // the density increases rapidly and both the direct and the 
conjugate Polyakov loops become nonzero, indicating "deconfinement" . 

The deviation from SU(3) during the complex Langevin evolution is shown 
in Fig. [151 using Trf/jf/4/3 as the observable. After the initial thermalization 
stage, this quantity fluctuates around 3.5 > 1. The fluctuations are similar for 
all values of the chemical potential we considered. Using spatial links f/j rather 
than 1/4^ gives a similar result. 
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Figure 14: Real part of the Polyakov loop (P) and the conjugate Polyakov loop 
(P^^) (left) and the density (n) (right) as a function of /x at /3 = 5.6, k = 0.12 
on a 4^ lattice, with Nf = 3 flavours. 




20000 30000 

Lansevin iteration 



Figure 15: Deviation from SU(3): TrUlU^/S during the Langevin evolution, for 
/i = 0.5 and 0.9. 



5.3 Phase of the determinant 

We study the phase of the determinant as before. Scatter plots during the 
Langevin evolution are shown in Fig. [161 for /^ = 0.5,0.6,0.7, and 0.8. At the 
smallest value of //, the average phase factor is close to one; for the real part we 
find 0.91 ± 0.28, while the imaginary part is consistent with zero (0.009 ± 0.39). 
At the larger values shown here, the distribution immediately becomes very wide 
and the average phase factor is consistent with zero (but with a large error). 
Note that the scale is very different compared to the one link model. Such an 
(apparently) abrupt change in the average phase factor when moving from a 
low-density to a high-density phase is somewhat reminiscent of what is found in 
random matrix studies, see e.g. Refs. [38 l 1391 1^ . 
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Figure 16: Scatter plot of e^*''^ = det M(/i)/det M(— /i) during the Langevin 
evolution for various values of /x at /3 = 5.6, k = 0.12 on a 4^ lattice. 

At large chemical potential the average phase factor approaches 1 again. This 
follows from the behaviour of the determinant, 



lim det M{fi) = Y\ /i^e^^/^ \l + 3e~''/^ (h + h'^) P^^ + 0(e"=^'^/^)l , (5.28) 

i^OO J. J. L 

X 

lim det M{-fi) = n h^e^^/^ \l + 'if^l^ (h + h^^) P, + 0{e^^^/^)] , (5.29) 

I — ^on XXL \ / J 



/^— >oo 



such that 

However, the values of the chemical potential we consider are not in that satura- 
tion region. 

6 Summary and outlook 

We have considered stochastic quantization for theories with a complex action 
due to finite chemical potential, and applied complex Langevin dynamics to U(l) 
and SU(3) one link models and QCD in the hopping expansion. In the latter, the 
full gauge dynamics is preserved but the fermion determinant is approximated. 
In all cases the complex determinant satisfies det M{fi) = [det M (— /i)]*, as is the 
case in QCD. We studied the (conjugate) Polyakov loops, the density and the 
phase of the determinant. In the one link models excellent agreement between 
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the numerical and exact results was obtained, for all values of /i ranging from 
zero to saturation. 

In the one link models the phase of the determinant was studied in detail. 
Even when the phase factor varies significantly during the Langevin evolution, 
its distribution is well-defined and expectation values can be evaluated without 
problems. The sign problem does not appear to be an obstruction here. In 
QCD in the hopping expansion, first results indicate that the behaviour of the 
average phase factor changes abruptly when moving from the low-density to a 
high-density phase. Nevertheless, other observables (Polyakov loop, density) are 
still under control, even at larger chemical potential. 

In the case of the U(l) model, we found strong hints why the sign problem does 
not appear to affect this method. We found that important features of classical 
flow and classical fixed points are largely independent of the chemical potential. 
The presence of /i changes the complex Langevin dynamics only quantitatively 
but not qualitatively, even when the average phase factor of the determinant 
becomes very small. Moreover, a study of the complex Fokker-Planck equation 
shows that all nonzero eigenvalues are real and positive, also in the presence 
of a nonzero chemical potential. An open question concerns the relationship 
between the stationary solution of the complex Fokker-Planck operator and its 
real counterpart. The structure in the U(l) model responsible for these results is 
related to symmetry properties of the action and the determinant. Therefore, it 
may be envisaged that they carry over to the more complicated cases. 

There are many directions into which this work can be extended. In the U(l) 
model, considerable insight could be obtained (semi-)analytically. It will be in- 
teresting to extend this analysis to more complicated theories. It will also be 
useful to perform further tests of the method in other simple models sharing rel- 
evant features with QCD at finite /i. Concerning QCD in the hopping expansion, 
for which we presented first results here, a more systematic study stays ahead. 
One way to test the approach is to also perform (real) Langevin dynamics at 
imaginary chemical potential, which goes smoothly and without runaway and 
convergence problems, and perform an analytical continuation. Finally, it will 
be interesting to apply this method to QCD both at large density as well as in 
the region of small chemical potential around the crossover temperature. Here 
it may shed light on the (non)existence of the critical point, in a setup which 
is manifestly independent from the other approaches available in this region. It 
should be noted that the Langevin method only requires the derivative of the 
determinant, and not the determinant itself. 
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A Fokker-Planck equation in Minkowski time 

In Refs. [26l [271 128] stochastic quantization was applied to nonequilibrium quan- 
tum dynamics in real time. For completeness, we give here the analysis of the 
corresponding complex Fokker-Planck equation for the one link U(l) modelo 
Consider the following partition function in Minkowski time, 

/'^ dx 
— e*^p, Sp = l3cosx+px. (A.l) 

-IT ^^ 

The term px, with p integer, is a reweighting term, used to stabilize the Langevin 
dynamics (see also Ref. [2H])- The Fokker-Planck equation for the (complex) 
distribution Pp reads 

d ^ , „^ d f d dS„\ „ , ^, 

-P,(,,B)=-[.--,^y,(.,B) 

= z/Pp (x, e)+i[p sin(x) - p] Pp(x, 9) + ip cos(x)Pp(x, 6). (A.2) 

Here ly is essentially the normalization of the noise: for z/ = 1 we have the full 
quantum case, for z/ = we obtain the classical evolution (which is equivalent to 
taking /3 and p to oo). 

To continue we discretize x a.s xi = 2nl/N, with — (A^ — l)/2 < / < A^/2, and 
define the modes as 

N/2 



Pp(n,e) = - ^ e-^'P,( 


xi,9), 


(A.S) 


l=-(N-l)/2 




N/2 




Pp{xu9)= )_J e-*"^'P,(n,^). 


(A.4) 


n=-(Af-l)/2 




The Fokker-Planck equation for the modes Pp(n, 9) then reads 




>,(»,») = - 


' 7V2 A^ " 
u siv?{kn/2)+p sm.{kn) 

TT^ //IT 


Pp{n,9) 











^This Appendix is partly based on Ref. [41]. 

33 



P N 

+1- — 

2 27r 



sin(/c„ 



+1 



)Pp{n 



sin(A;. 



n~l 



)Pp(n-l, 



./5 



PJn + 1, 1 



PM 



where kn = 2'Kn/N. For small n/N this reduces to 



^/pin, 



— [vn 



+ pn)Pp{n,i 



2 ^ 



Pp(n 



Pp(n 



(A.5) 



(A.6) 



which can be obtained directly from the continuum Fokker-Planck equation before 
discretizing x. Extension to general (complex) (3 = (3r + i[3i and p = Pr + ipi is 
straightforward. In the following no explicit u means u = 1. 
From averages with the distribution P we obtain 



{e'n^ 



{e'n^ 



f:jxe^^^Pp{x,9) _ P^(q,9) 

j:jxP,ix,e) p,(o,^)' 

/^^ dx e^^^Pola;, 9) Z^ dx e^('?-P)^Pp(x, 9) 



r_JxPo( 



X, 



Zjxe-'P-Ppix,9) 



Pp{q-P,0] 
Pp{-P,0) 



(A.7) 
(A.8) 



Notice that p and q (both integers) are interchangeable. This implies that simu- 
lations can be performed at p = 0, while Eqs. (1A.7I) - ( 1A.8I) can be used to obtain 
Pp for any p. 

Similar to the procedure of Sec. 13. 4[ we have solved the complex Fokker- 
Planck equation numerically for the modes Pp{n,9). In contrast to the finite /x 
case, these modes are now complex in general. For example, for even (odd) p 
even (odd) modes are real and odd (even) modes are imaginary, in agreement 
with the symmetries of the action. We further find that for p > the solution for 
positive modes converges correctly, but negative modes diverge, and vice versa. 
The numerical solution for p = converges quickly to the values determined by 
e*'^°, when < /9 < 2.3|j In Fig. [17] (left) we show the Langevin time dependence 
of some modes when (3 = 1, p = 0. Using the asymptotic values for the modes 
and Eqs. (1A.7P - (lA.Sp . we can reconstruct Pp(x), which agrees nicely with e''^^'-'^). 
The Langevin simulation itself also yields good results for the expectation values 
when p 7^ or at p = 0, provided reweighting from p 7^ is used. For a more 
thorough discussion of the conditions for convergence of the Langevin simulation, 
see Ref. [2H]. 

Again further insight can be obtained from the eigenvalues, determined by 
the eigenvalue equation 



nin 



p)Pp{n) 



n 



-ip 



Pp{n 



Pp{n 



\Pp{n). 



(A.9) 



If A 7^ 0, Pp(0) vanishes and the sequences for n > and n < split again in 
two. Positive and negative n are related by changing the sign of p. In Fig. [T7| we 



'Notice that the partition function Zo has a zero at (3 near 2.4. 
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Figure 17: Complex Fokker-Planck equation for Minkowski dynamics. Left: 
Langevin time dependence of the modes Po{n,6) for various values of n, with 
(3 = 1, p = 0. Right: Smallest nonzero eigenvalue of the complex Fokker-Planck 
equation as a function of /3 for various values of the reweighting parameter p. 

show the smallest nonzero eigenvalue for the positive n sequence for five values 
of p. In contrast to the finite /i case, we find that eigenvalues may be negative, 
depending on the value of p and (3, indicating the possibility of problems with 
convergence and stability. This corroborates the above observations. 

Finally, one may also study the real Fokker-Planck equation to obtain the 
true probability distribution pp{x,y,9) and its modes. For an analytic function 
0{z) we have 



dxdy 



27r 



Pp{x,y,9)0{x + iy) 



hence, in particular 

dxdy 



2n 



X, y, ^)e*"--"^ 



dye 



""^Pp(n, y, 



(IT 

^Pp{x,e)0{x), 



^PA^.O)e^ 



(A.IO) 



Pp{n, 



(A.ll) 
The expectation values with pp should represent the averages over the Langevin 
process itself. Even when the latter converge to the correct values, the real 
Fokker-Planck equation does not show good behaviour, however. We thus have 
the situation that we find agreement between the complex Fokker-Planck equation 
(with the corresponding complex distribution P) and the actual Langevin process, 
while the true probability distribution pp, which formally mediates between the 
former two, appears more difficult to control. 
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